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ABSTRACT 


Numerical  Methods  are  developed  for  constructing  body-fitted 
curvilinear  coordinate  system  for  the  region  surrounding  an  arbitrary 
t hree-d laens lonal  body.  Finite  difference  acheaes  are  Investigated 
for  solving  the  Havler-Stokes  equations  on  body-fitted  coordinates. 
Solutions  are  computed  for  flow  about  a  sphere  and  a  finite  wing. 
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1.  Introduction 


The  research  efforts  of  the  past  three  years  on  three-dimensional 
laminar  flow  can  be  divided  into  two  areas.  The  first  task  was  the 
generation  of  a  body-fitted  coordinate  system  by  solving  a  system  of 
elliptic  equations.  The  theoretical  basis  for  the  coordinate  generation 
scheme  appeared  in  a  paper  by  Mastln  and  Thompson  [1].  Although  the 
method  was  developed  for  the  region  about  a  single  body,  the  capability 
exists  for  generating  coordinate  systems  for  nearly  any  bounded  three- 
dimensional  region.  It  may  be  necessary  to  partition  the  physical 
region  into  simpler  subregions.  This  partitioning,  which  is  essential 
for  very  complicated  regions,  has  proven  to  be  easily  Implemented  in 
the  coordinate  generation  scheme. 

The  second  area  of  research,  which  depended  on  the  generation  of 
a  suitable  coordinate  system,  was  the  numerical  solution  of  the  Navier- 
Stokes  equations  for  viscous  flow  about  an  arbitrary  body.  Two  Implicit 
methods  have  been  developed.  Both  are  generalizations  to  curvilinear 
coordinates  of  methods  which  have  been  developed  for  rectangular  coordinate 
systems.  The  use  of  SOR  iteration  in  solving  the  Navier-Stokes  equations 
on  three-dimensional  body-fitted  coordinates  was  first  reported  by  Mastln 
and  Thompson  [2].  More  extensive  computations  and  comparisons  with  known 
experimental  results  for  viscous  flow  about  a  sphere  will  appear  in  the 
dissertation  by  Fu  [3]. 

In  an  effort  to  develop  a  more  efficient  method  on  the  larger  compu¬ 
tational  fields  needed  for  high  Reynolds  number  calculations,  a  fractional 


step  aethod  was  Investigated.  The  basic  algorithm  and  preliminary 
computational  results  have  appeared  In  a  paper  by  Mastin,  Ghosh,  and 
Thompson  [4].  In  the  problems  considered  so  far,  the  physical  region 
was  partitioned  and  hence,  an  Iteration  process  was  needed  with  the 
fractional  step  algorithm  so  that  the  equations  are  satisfied  across 
the  cuts.  While  this  allows  the  user  to  handle  larger  and  more  complex 
regions,  it  decreased  the  efficiency  of  the  method.  The  accuracy  of 
the  fractional  step  method  has  been  about  the  same  as  the  SOR  method 
for  flow  about  a  sphere. 

The  following  sections  will  describe  the  major  accomplishments 
of  this  research  project.  A  more  detailed  description  can  be  found 
by  referring  to  the  references  given  above. 

2.  Body-Fitted  Coordinate  System 

Curvilinear  coordinate  systems  have  been  generated  about  various 
three-dimensional  bodies.  The  physical  region  is  truncated  at  a  finite 
distance  and  the  region  between  the  surface  of  the  body  and  the  surface 
at  infinity  is  partitioned  into  subregions  which  can  be  transformed  to 
rectangular  computational  regions.  For  low  Reynolds  number  problems, 
three  subregions  were  used  as  in  Figure  1.  For  higher  Reynolds  numbers, 
the  central  section  was  divided  into  four  smaller  subregions  to  give  a 
total  of  six  subregions  of  approximately  the  same  size.  The  partitioning 
of  the  region  served  a  dual  purpose.  First  of  all,  it  allowed  for  a 
nearly  equal  distribution  of  mesh  points  on  the  surface  of  the  body. 

This  would  not  be  the  case  for  a  spherical  or  cylindrical  type  coordinate 
system  where  mesh  points  would  be  clustered  along  an  axis  where  the 
transformation  would  be  singular.  Secondly,  it  allowed  for  a  much  larger 
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computational  mesh  since  it  was  only  necessary  to  keep  the  data  for  one 
region  in  the  core  memory  of  the  computer. 

The  transformation  from  the  physical  xyz-space  to  the  computational 
yvC-apace  was  obtained  by  solving  the  elliptic  system 

72p  ■  f(ii,v.O 

?2v  -  g(u,v,5) 

V25  ■  h(y,v,£)  (1) 

The  homogeneous  equations  (f  ■  g  *  h  ■  0)  give  a  smooth  transformation  with 

a  nonvanishing  Jacobian.  In  some  cases  a  nonzero  value  of  h  was  used  so 
that  more  mesh  points  were  concentrated  near  the  body.  The  generality 
of  the  transformation  method  has  been  verified  by  constructing  coordinate 
systems  about  a  multitude  of  bodies  of  varying  shapes.  These  Include  a 
sphere,  ellipsoids,  cylinders  with  hemispherical  caps,  cylinders  with 
conical  caps,  finite  wings,  and  a  flat  plate. 

A  difficult  problem  in  the  generation  of  coordinate  systems  for  three- 
dimensional  bodies  is  the  specification  of  the  boundary  correspondence  for 
the  system  (1).  In  order  to  simplify  the  determination  of  coordinates  for 
mash  points  on  the  body,  these  mesh  points  were  chosen  to  lie  on  cross- 
sections  of  the  body.  This  also  aided  in  the  plotting  of  the  flow  field 
data. 

The  initial  velocity  and  pressure  distribution  for  the  solution  of  the 
Navier-Stokes  equations  were  usually  obtained  from  a  potential  flow  calcu¬ 
lation.  The  potential  function  can  be  easily  calculated  during  the  solution 
of  the  system  (1)  in  the  computational  regions. 
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thods  for  Cho  Navler-Stokes  Equations 


Successive  over relexet Ion  (SOR)  iterative  Methods  have  been  successfully 
used  for  solving  the  Navler-Stokes  equations  on  nany  two-dinsnalonal  regions 
for  a  wide  range  of  Reynolds  numbers.  A  straightforward  generalisation  to 
three  dlaensions  has  worked  equally  well  although  the  Method  has  proven 
expensive  in  terns  of  conputer  tine  and  storage. 

The  form  of  the  Navler-Stokes  equations  most  frequently  used  was 

V  •  v  -  0  (2a) 

v  +  (v  •  V)v  ■  -  Vp  +  ^  V2v  (2b) 

where  v  •  ut  +  vj  +  w£  is  the  velocity,  p  is  the  pressure,  and  R  is  the 
Reynolds  nunber.  A  few  calculations  were  nade  with  the  advectlve  terms  in 
conservation  form.  The  result  was  a  longer  run  time  with  no  appreciable 
difference  in  the  solution.  The  most  efficient  manner  of  differencing  the 
momentum  equation  (2b)  was 

r*1  -  ^-VtK^.V)?”  +  Vp”*1  -  i  V*?”*1] .  (3) 

where  n  is  the  tine  step  index  and  Vt  the  step  size.  In  this  formulation, 
the  viscous  and  pressure  tens  are  Implicit  while  the  advectlve  terns  are 
explicit.  The  stability  criteria  imposed  by  the  explicit  terms  was  quite 
mild  due  to  the  snail  velocity  components  near  the  body  where  the  mesh 


spacing  was  also  snail.  In  fact,  with  a  fully  lnpllcit  scheme  which  was 
also  coded,  a  smaller  time  step  size  was  needed  to  maintain  SOR  convergence 
than  was  required  for  the  stability  of  the  above  scheme. 

The  Poisson  equation  used  for  the  pressure  calculation  was 

'V1  -  <»v  <«> 
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where  D  -  V  •  v.  It  «u  necessary  to  delate  sobs  terms  In  the  pressure 
equation  derived  from  (3)  so  that  the  resulting  systen  of  four  equations 
could  be  solved  by  SOR  iteration  in  the  transformed  regions  to  obtain  values 
of  v  end  p  at  tine  step  n+1. 

In  ell  computations  a  unit  free  stream  velocity  was  imposed  on  the 
surface  at  infinity.  On  the  body  a  no-slip  velocity  condition  was  assumed 
and  a  normal  derivative  condition  for  the  pressure  was  computed  from  the 
equation 

Vp-*1  -  i  7^1  (5) 

A  constant  pressure  on  the  surface  at  infinity  was  used  for  most  of  the 
computations.  One  could  argue  that  the  same  normal  derivative  condition 
should  hold  on  the  surface  at  infinity.  Thus  condition  was  tried  and  gave 
essentially  the  same  pressure  distribution  on  the  body.  It  does,  however, 
require  the  determination  of  the  pressure  level  at  Infinity  in  order  to 
plot  the  pressure  coefficient. 

The  above  schema  gave  realistic  solutions  to  fluid  flow  problems  and 
the  velocity  and  pressure  values  compared  well  with  known  analytic  and 
nimMrlcal  results.  The  calculated  value  of  D  was  large  near  the  body  where 
an  impulsive  start  was  used  and  decreased  slowly  as  time  increased.  It 
was  also  observed  that  the  velocity  values  converged  much  faster  than  the 
pressure  when  the  system  was  solved  by  SOR.  Consequently,  a  second  SOR 
method  was  developed  to  solve  the  Navier-Stokes  equations.  An  auxiliary 

-jxi /2 

vector  v  '  was  introduced  in  two  steps. 


-m+1/2  -*n  r  ,-*n  _.-ni  1  _2-*n+l/2. 

v  -  v  -  Vt[(v  •  V)v  -  V*v  ] 


V*1  -  ^+1/2  -  vt  Vp1*1 


(6a) 

(6b) 


- - - - - - - - -  - - 


The  pressure  is  computed  sffcsr  the  first  step  by  solving  the  Poisson 
equation 

.  JL  ,  .  r*in  .  (7) 

As  with  many  splitting  techniques ,  Instabilities  esn  arise  if  correct 
boundary  values  are  not  chosen  for  the  auxiliary  functions  computed  at  the 
intermediate  steps.  For  simplicity  In  programming,  the  first  value  used  was 

£*+1/2  .  £n 

on  the  boundary.  This  value  worked  well  In  all  our  calculations  with  no 
Indication  of  Instability.  Two  normal  derivative  conditions  for  pressure 
were  used.  Both  were  obtained  by  considering  the  limiting  value  of  the 
pressure  gradient  at  the  surface  of  the  body.  If  only  the  second  equation 
(6b)  is  considered  we  have 

Vp®*1  -  0  ,  (8) 

while  if  both  equations  are  considered  we  arlve  at 

vPM1 .  i  .  <9> 

For  this  type  of  splitting,  the  first  alternative  (8)  has  been  used  by 
others  in  the  theoretical  analysis  of  the  method  and  has  also  been  used  in 
most  of  our  calculations.  The  normal  component  of  the  pressure  gradient 
in  (9)  would  appear  to  give  a  solution  in  closer  agreement  with  the  unspllt 
method  where  (5)  is  used.  However,  both  boundary  conditions  gave  essentially 
the  same  results  in  the  problems  which  have  been  considered  in  our  work. 

This  two  step  method  was  successfully  applied  and  when  the  results 
of  both  methods  were  compared,  the  difference  between  the  two  solutions 
was  less  than  the  overall  error  in  the  approximation  of  the  true  solution. 

The  two  step  method  was  preferred  because  it  was  faster  and  gave  a  smaller 
value  of  the  divergence  D.  Unfortunately,  neither  method  was  able  to 


1 


reduce  Che  maximum  absolute  value  of  D  below  about  10. Vt  for  aost  probleas. 
For  high  Reynolds  nuaber  calculations,  the  two  step  method  was  used  with  a 
sero  noraal  pressure  derivative.  The  one  step  net hod  did  not  work  as  well 
on  the  highly  contracted  coordinate  systeas  needed  for  high  Reynolds  nuaber 
calculations. 


4.  A  Fractional  Step  Method  for  the  Navler-Stokes  Equations 

The  fractional  step  nethod  used  in  this  project  involved  a  splitting 
of  the  pressure  tern  from  the  moaentua  equations  as  in  the  two  step  method 
discussed  earlier  and  a  splitting  of  the  reaalning  terms  according  to  the 
derivative  direction  in  computational  space.  In  order  to  solve  the  pressure 
equation  in  the  same  manner,  a  parabolic  pressure  equation  is  needed.  Thus 
the  continuity  equation  (2a)  is  replaced  by  the  artificial  compressibility 


equation. 


ept  +  V  •  v  "  0 


It  was  discovered  that  better  results  could  be  obtained  by  taking  a 
variable  e  which  was  0(At2). 

The  same  type  of  partitioning  was  used  in  the  physical  region  as  with 
the  SOR  methods.  An  iteration  scheme  was  Introduced  at  each  fractional  step 
so  that  all  equations  are  satisfied  on  the  surfaces  which  partition  the 
region.  In  this  Iteration  process  the  nonlinear  coefficients  and  mixed 
derivative  terms  are  also  updated  so  that  all  terms  are  treated  implicitly. 

The  splitting  of  the  velocity  derivative  terms  in  the  momentum  equations 
according  to  coordinate  directions  in  the  computational  regions  can  be 


expressed  as 


I  o-*  4  4 

(v  .  7)v  -  ■=■  V2v  ■  A  v  +  A  v  +  A_v  . 

R  p  v  C 


7 


e(pn+2/3  _  pO+1/3)  .  (At)2  *y*2/3  -  0  (13b) 

etp"*-1  -  p“+2/3)  +  At (V.v) 0+3/4  -  (At)2  B^pn+1  -  0  (13c) 

The  divergence  term  has  been  included  in  the  last  equation  for  simplicity. 

As  with  the  two  step  SOR  scheme,  there  are  two  possible  boundary  conditions 
for  the  body  pressure  which  can  be  derived  by  considering  (12d)  alone  or 
by  summing  (12a)-(12d).  The  above  equations  (13a)- (13c)  are  solved  by  a 
tridiagonal  algorithm  in  an  iterative  sequence  which  is  used  to  assure  that 
all  equations  are  satisfied  on  the  surfaces  which  partition  the  region. 

The  pressure  boundary  condition  is  also  applied  in  the  iterative  solution 
of  (13c).  Even  for  a  single  computational  region  with  no  partitioning  it 
may  be  necessary  to  implement  the  normal  derivative  condition  for  the 
pressure  function  iteratively  to  avoid  stability  problems. 

This  scheme  has  been  used  to  compute  viscous  flow  about  a  sphere  at 
low  and  intermediate  Reynolds  numbers.  It  is  comparable  to  the  SOR 
methods  in  both  accuracy  and  computer  time  used.  However,  the  full 
potential  of  the  method  has  not  been  realized  in  the  present  configuration 
with  the  partitioned  regions.  A  direct  Implementation  of  the  method, 
with  no  iterations,  would  be  an  order  of  magnitude  faster.  The  fractional 
step  method  also  has  an  advantage  when  core  storage  is  considered.  At 
each  fractional  step  only  one  third  of  the  derivative  coefficients  of  the 
transformed  equations  are  required  for  the  computations.  The  main  difficulty 
baa  been  the  choice  of  a  proper  value  for  the  parameter  e.  A  value  of  c 
which  is  too  large  leaves  the  pressure  practically  unchanged  and  the 
velocity  divergence  becomes  large.  On  the  other  hand,  a  value  of  e  which 
is  too  small  causes  the  solution  to  become  unbounded  after  very  few  time 
steps. 
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utatlooal  Results 


The  principle  objective  of  thia  reaearch  has  been  algorithm  development 


rather  than  the  developaent  of  production  codes.  Thus  nost  of  the  Navier 


Stokes  solutions  have  been  computed  for  flow  about  a  sphere  where  valid 


experimental  and  numerical  results  are  available  for  comparison.  The 


present  methods  worked  best  for  the  Intermediate  Reynolds  numbers  between 


40  and  1,000,  however,  solutions  were  obtained  for  Reynolds  numbers 


The  coordinate  system  about  a  sphere  was  constructed  using  the  schei 


described  in  Section  2.  Figures  2  and  3  contain  plots  of  a  crossectlon  of 


the  coordinte  surfaces  snd  a  perspective  plot  of  one  coordinate  surface 


surrounding  the  sphere.  For  viscous  flow  about  a  sphere,  drag  calculations 
compared  well  with  experimental  results  as  can  be  seen  in  Figure  4.  Both 


the  SOR  and  fractional  step  schemes  gave  clearly  formed  vortices  at  the 


rear  of  the  sphere.  A  crossectlon  of  the  steady-state  velocity  field  for 


the  SOR  computation  at  a  Reynolds  number  of  40  appears  in  Figure  5.  The 


corresponding  pressure  coefficient  is  plotted  in  Figure  6.  As  the  Reynolds 


number  increases,  the  velocity  in  the  vortex  region  Increases.  This  can  be 


seen  by  examining  Figure  7  which  shows  the  same  crossectlon  with  a  Reynolds 


number  of  290.  The  vectors  have  been  scaled  to  one- fourth  their  actual 


length.  The  drop  in  pressure  at  the  forward  stagnation  point  can  be  noted 


from  Figure  8.  This  second  solution  was  computed  using  the  fractional  step 


The  practical  application  of  body-fitted  coordinate  systems  has  not 
been  ignored  in  our  research.  Flow  about  finite  wings  have  been  computed 


for  wings  having  various  lengths,  cord  to  thickness  ratios,  camber,  and 


tapering.  Some  of  the  large  scale  characteristics  of  the  flow  can  be 


observed  even  though  the  wake  region  was  not  properly  modeled  due  to  the 
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relatively  small  number  of  mesh  points  downstream  of  the  wing.  For  example. 


consider  the  wing  in  Figure  9  which  is  perpendicular  to  the  free  stream  and 
at  a  10  degree  angle  of  attack.  An  SOR  scheme  was  used  with  a  Reynolds 
number  of  40.  A  rear  view  of  the  velocity  near  the  trailing  edge,  illustrated 


In  Figure  10,  indicates  a  vortex  formed  at  the  wing  tip  and  downwash 
immediately  behind  the  wing.  Due  to  symmetry,  only  half  the  wing  was 


plotted. 


6.  Conclusions 


Body-fitted  coordinate  systems  have  proven  to  be  a  valuable  tool  for 


solving  the  Navier-Stokes  equations  for  viscons  flow  about  arbitrary  bodies. 
Both  SOR  and  fractional  step  methods  have  been  implemented  on  body-fitted 


coordinates.  Thus  the  programs  for  solving  the  Navier-Stokes  equations  are 


essentially  Independent  of  the  coordinate  system  even  if  several  rectangular 
computational  regions  are  used.  The  accuracy  of  the  methods  has  been  tested 
by  computing  viscous  flow  about  a  sphere  for  a  wide  range  of  Reynolds  numbers . 


Despite  the  ability  to  use  relatively  large  time  steps,  the  solution  of  the 


Navier-Stokes  equations  starting  with  potential  flew  and  marching  to  the 


steady-state  solution  has  required  large  amounts  of  computer  time.  However, 


realistic  solutions  of  viscous  flow  problems  about  a  variety  of  body  shapes 


and  a  wide  range  of  Reynolds  numbers  can  be  computed  by  the  methods  developed 


under  this  project. 
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Figure  4.  Drag  coefficient  for  spheres  as  a  function  of  the  Reynolds  number 
Curve  (1):  Stokes's  theory,  curve  (2):  Oseen’a  theory. 
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